home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / fortran / linpklib.zip / SSVDC.FOR < prev    next >
Text File  |  1984-01-06  |  15KB  |  481 lines

  1.       SUBROUTINE SSVDC(X,LDX,N,P,S,E,U,LDU,V,LDV,WORK,JOB,INFO)
  2.       INTEGER LDX,N,P,LDU,LDV,JOB,INFO
  3.       REAL X(LDX,1),S(1),E(1),U(LDU,1),V(LDV,1),WORK(1)
  4. C
  5. C
  6. C     SSVDC IS A SUBROUTINE TO REDUCE A REAL NXP MATRIX X BY
  7. C     ORTHOGONAL TRANSFORMATIONS U AND V TO DIAGONAL FORM.  THE
  8. C     DIAGONAL ELEMENTS S(I) ARE THE SINGULAR VALUES OF X.  THE
  9. C     COLUMNS OF U ARE THE CORRESPONDING LEFT SINGULAR VECTORS,
  10. C     AND THE COLUMNS OF V THE RIGHT SINGULAR VECTORS.
  11. C
  12. C     ON ENTRY
  13. C
  14. C         X         REAL(LDX,P), WHERE LDX.GE.N.
  15. C                   X CONTAINS THE MATRIX WHOSE SINGULAR VALUE
  16. C                   DECOMPOSITION IS TO BE COMPUTED.  X IS
  17. C                   DESTROYED BY SSVDC.
  18. C
  19. C         LDX       INTEGER.
  20. C                   LDX IS THE LEADING DIMENSION OF THE ARRAY X.
  21. C
  22. C         N         INTEGER.
  23. C                   N IS THE NUMBER OF COLUMNS OF THE MATRIX X.
  24. C
  25. C         P         INTEGER.
  26. C                   P IS THE NUMBER OF ROWS OF THE MATRIX X.
  27. C
  28. C         LDU       INTEGER.
  29. C                   LDU IS THE LEADING DIMENSION OF THE ARRAY U.
  30. C                   (SEE BELOW).
  31. C
  32. C         LDV       INTEGER.
  33. C                   LDV IS THE LEADING DIMENSION OF THE ARRAY V.
  34. C                   (SEE BELOW).
  35. C
  36. C         WORK      REAL(N).
  37. C                   WORK IS A SCRATCH ARRAY.
  38. C
  39. C         JOB       INTEGER.
  40. C                   JOB CONTROLS THE COMPUTATION OF THE SINGULAR
  41. C                   VECTORS.  IT HAS THE DECIMAL EXPANSION AB
  42. C                   WITH THE FOLLOWING MEANING
  43. C
  44. C                        A.EQ.0    DO NOT COMPUTE THE LEFT SINGULAR
  45. C                                  VECTORS.
  46. C                        A.EQ.1    RETURN THE N LEFT SINGULAR VECTORS
  47. C                                  IN U.
  48. C                        A.GE.2    RETURN THE FIRST MIN(N,P) SINGULAR
  49. C                                  VECTORS IN U.
  50. C                        B.EQ.0    DO NOT COMPUTE THE RIGHT SINGULAR
  51. C                                  VECTORS.
  52. C                        B.EQ.1    RETURN THE RIGHT SINGULAR VECTORS
  53. C                                  IN V.
  54. C
  55. C     ON RETURN
  56. C
  57. C         S         REAL(MM), WHERE MM=MIN(N+1,P).
  58. C                   THE FIRST MIN(N,P) ENTRIES OF S CONTAIN THE
  59. C                   SINGULAR VALUES OF X ARRANGED IN DESCENDING
  60. C                   ORDER OF MAGNITUDE.
  61. C
  62. C         E         REAL(P).
  63. C                   E ORDINARILY CONTAINS ZEROS.  HOWEVER SEE THE
  64. C                   DISCUSSION OF INFO FOR EXCEPTIONS.
  65. C
  66. C         U         REAL(LDU,K), WHERE LDU.GE.N.  IF JOBA.EQ.1 THEN
  67. C                                   K.EQ.N, IF JOBA.GE.2 THEN
  68. C                                   K.EQ.MIN(N,P).
  69. C                   U CONTAINS THE MATRIX OF RIGHT SINGULAR VECTORS.
  70. C                   U IS NOT REFERENCED IF JOBA.EQ.0.  IF N.LE.P
  71. C                   OR IF JOBA.EQ.2, THEN U MAY BE IDENTIFIED WITH X
  72. C                   IN THE SUBROUTINE CALL.
  73. C
  74. C         V         REAL(LDV,P), WHERE LDV.GE.P.
  75. C                   V CONTAINS THE MATRIX OF RIGHT SINGULAR VECTORS.
  76. C                   V IS NOT REFERENCED IF JOB.EQ.0.  IF P.LE.N,
  77. C                   THEN V MAY BE IDENTIFIED WITH X IN THE
  78. C                   SUBROUTINE CALL.
  79. C
  80. C         INFO      INTEGER.
  81. C                   THE SINGULAR VALUES (AND THEIR CORRESPONDING
  82. C                   SINGULAR VECTORS) S(INFO+1),S(INFO+2),...,S(M)
  83. C                   ARE CORRECT (HERE M=MIN(N,P)).  THUS IF
  84. C                   INFO.EQ.0, ALL THE SINGULAR VALUES AND THEIR
  85. C                   VECTORS ARE CORRECT.  IN ANY EVENT, THE MATRIX
  86. C                   B = TRANS(U)*X*V IS THE BIDIAGONAL MATRIX
  87. C                   WITH THE ELEMENTS OF S ON ITS DIAGONAL AND THE
  88. C                   ELEMENTS OF E ON ITS SUPER-DIAGONAL (TRANS(U)
  89. C                   IS THE TRANSPOSE OF U).  THUS THE SINGULAR
  90. C                   VALUES OF X AND B ARE THE SAME.
  91. C
  92. C     LINPACK. THIS VERSION DATED 03/19/79 .
  93. C     G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
  94. C
  95. C     ***** USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS.
  96. C
  97. C     EXTERNAL SROT
  98. C     BLAS SAXPY,SDOT,SSCAL,SSWAP,SNRM2,SROTG
  99. C     FORTRAN ABS,AMAX1,MAX0,MIN0,MOD,SQRT
  100. C
  101. C     INTERNAL VARIABLES
  102. C
  103.       INTEGER I,ITER,J,JOBU,K,KASE,KK,L,LL,LLS,LM1,LP1,LS,LU,M,MAXIT,
  104.      *        MM,MM1,MP1,NCT,NCTP1,NCU,NRT,NRTP1
  105.       REAL SDOT,T,R
  106.       REAL B,C,CS,EL,EMM1,F,G,SNRM2,SCALE,SHIFT,SL,SM,SN,SMM1,T1,TEST,
  107.      *     ZTEST
  108.       LOGICAL WANTU,WANTV
  109. C
  110. C
  111. C     SET THE MAXIMUM NUMBER OF ITERATIONS.
  112. C
  113.       MAXIT = 30
  114. C
  115. C     DETERMINE WHAT IS TO BE COMPUTED.
  116. C
  117.       WANTU = .FALSE.
  118.       WANTV = .FALSE.
  119.       JOBU = MOD(JOB,100)/10
  120.       NCU = N
  121.       IF (JOBU .GT. 1) NCU = MIN0(N,P)
  122.       IF (JOBU .NE. 0) WANTU = .TRUE.
  123.       IF (MOD(JOB,10) .NE. 0) WANTV = .TRUE.
  124. C
  125. C     REDUCE X TO BIDIAGONAL FORM, STORING THE DIAGONAL ELEMENTS
  126. C     IN S AND THE SUPER-DIAGONAL ELEMENTS IN E.
  127. C
  128.       INFO = 0
  129.       NCT = MIN0(N-1,P)
  130.       NRT = MAX0(0,MIN0(P-2,N))
  131.       LU = MAX0(NCT,NRT)
  132.       IF (LU .LT. 1) GO TO 170
  133.       DO 160 L = 1, LU
  134.          LP1 = L + 1
  135.          IF (L .GT. NCT) GO TO 20
  136. C
  137. C           COMPUTE THE TRANSFORMATION FOR THE L-TH COLUMN AND
  138. C           PLACE THE L-TH DIAGONAL IN S(L).
  139. C
  140.             S(L) = SNRM2(N-L+1,X(L,L),1)
  141.             IF (S(L) .EQ. 0.0E0) GO TO 10
  142.                IF (X(L,L) .NE. 0.0E0) S(L) = SIGN(S(L),X(L,L))
  143.                CALL SSCAL(N-L+1,1.0E0/S(L),X(L,L),1)
  144.                X(L,L) = 1.0E0 + X(L,L)
  145.    10       CONTINUE
  146.             S(L) = -S(L)
  147.    20    CONTINUE
  148.          IF (P .LT. LP1) GO TO 50
  149.          DO 40 J = LP1, P
  150.             IF (L .GT. NCT) GO TO 30
  151.             IF (S(L) .EQ. 0.0E0) GO TO 30
  152. C
  153. C              APPLY THE TRANSFORMATION.
  154. C
  155.                T = -SDOT(N-L+1,X(L,L),1,X(L,J),1)/X(L,L)
  156.                CALL SAXPY(N-L+1,T,X(L,L),1,X(L,J),1)
  157.    30       CONTINUE
  158. C
  159. C           PLACE THE L-TH ROW OF X INTO  E FOR THE
  160. C           SUBSEQUENT CALCULATION OF THE ROW TRANSFORMATION.
  161. C
  162.             E(J) = X(L,J)
  163.    40    CONTINUE
  164.    50    CONTINUE
  165.          IF (.NOT.WANTU .OR. L .GT. NCT) GO TO 70
  166. C
  167. C           PLACE THE TRANSFORMATION IN U FOR SUBSEQUENT BACK
  168. C           MULTIPLICATION.
  169. C
  170.             DO 60 I = L, N
  171.                U(I,L) = X(I,L)
  172.    60       CONTINUE
  173.    70    CONTINUE
  174.          IF (L .GT. NRT) GO TO 150
  175. C
  176. C           COMPUTE THE L-TH ROW TRANSFORMATION AND PLACE THE
  177. C           L-TH SUPER-DIAGONAL IN E(L).
  178. C
  179.             E(L) = SNRM2(P-L,E(LP1),1)
  180.             IF (E(L) .EQ. 0.0E0) GO TO 80
  181.                IF (E(LP1) .NE. 0.0E0) E(L) = SIGN(E(L),E(LP1))
  182.                CALL SSCAL(P-L,1.0E0/E(L),E(LP1),1)
  183.                E(LP1) = 1.0E0 + E(LP1)
  184.    80       CONTINUE
  185.             E(L) = -E(L)
  186.             IF (LP1 .GT. N .OR. E(L) .EQ. 0.0E0) GO TO 120
  187. C
  188. C              APPLY THE TRANSFORMATION.
  189. C
  190.                DO 90 I = LP1, N
  191.                   WORK(I) = 0.0E0
  192.    90          CONTINUE
  193.                DO 100 J = LP1, P
  194.                   CALL SAXPY(N-L,E(J),X(LP1,J),1,WORK(LP1),1)
  195.   100          CONTINUE
  196.                DO 110 J = LP1, P
  197.                   CALL SAXPY(N-L,-E(J)/E(LP1),WORK(LP1),1,X(LP1,J),1)
  198.   110          CONTINUE
  199.   120       CONTINUE
  200.             IF (.NOT.WANTV) GO TO 140
  201. C
  202. C              PLACE THE TRANSFORMATION IN V FOR SUBSEQUENT
  203. C              BACK MULTIPLICATION.
  204. C
  205.                DO 130 I = LP1, P
  206.                   V(I,L) = E(I)
  207.   130          CONTINUE
  208.   140       CONTINUE
  209.   150    CONTINUE
  210.   160 CONTINUE
  211.   170 CONTINUE
  212. C
  213. C     SET UP THE FINAL BIDIAGONAL MATRIX OR ORDER M.
  214. C
  215.       M = MIN0(P,N+1)
  216.       NCTP1 = NCT + 1
  217.       NRTP1 = NRT + 1
  218.       IF (NCT .LT. P) S(NCTP1) = X(NCTP1,NCTP1)
  219.       IF (N .LT. M) S(M) = 0.0E0
  220.       IF (NRTP1 .LT. M) E(NRTP1) = X(NRTP1,M)
  221.       E(M) = 0.0E0
  222. C
  223. C     IF REQUIRED, GENERATE U.
  224. C
  225.       IF (.NOT.WANTU) GO TO 300
  226.          IF (NCU .LT. NCTP1) GO TO 200
  227.          DO 190 J = NCTP1, NCU
  228.             DO 180 I = 1, N
  229.                U(I,J) = 0.0E0
  230.   180       CONTINUE
  231.             U(J,J) = 1.0E0
  232.   190    CONTINUE
  233.   200    CONTINUE
  234.          IF (NCT